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Abstract. Climate models have traditionally only repre- 
sented heat and water fluxes within relatively shallow soil 
layers, but there is increasing interest in the possible role of 
heat and water exchanges with the deeper subsurface. Here, 
we integrate an idealized 50 m deep aquifer into the land sur- 
face module of the GISS ModelE general circulation model 
to test the influence of aquifer-soil moisture and heat ex- 
changes on climate variables. We evaluate the impact on the 
modeled climate of aquifer-soil heat and water fluxes sepa- 
rately, as well as in combination. The addition of the aquifer 
to ModelE has limited impact on annual-mean climate, with 
little change in global mean land temperature, precipitation, 
or evaporation. The seasonal amplitude of deep soil tempera- 
ture is strongly damped by the soil-aquifer heat flux. This not 
only improves the model representation of permafrost area 
but propagates to the surface, resulting in an increase in the 
seasonal amplitude of surface air temperature of > 1 K in the 
Arctic. The soil-aquifer water and heat fluxes both slightly 
decrease interannual variability in soil moisture and in land- 
surface temperature, and decrease the soil moisture memory 
of the land surface on seasonal to annual timescales. The re- 
sults of this experiment suggest that deepening the modeled 
land surface, compared to modeling only a shallower soil col- 
umn with a no-flux bottom boundary condition, has limited 
impact on mean climate but does affect seasonality and inter- 
annual persistence. 


1 Introduction 


The land- surface components of global climate models 
(GCMs) have typically represented water and heat storage 
as occurring only in a surface layer 1-5 m deep (Manabe, 
1969; Koster and Suarez, 1996; Milly and Shmakin, 2002; 
Essery et al., 2003; Dickinson et ah, 2006; Schmidt et ah, 
2006). However, deeper aquifers that contain groundwater 
form a critical and depleting portion of the exploitable water 
resource in many areas (Siebert et ah, 2010). These aquifers 
can exchange moisture with shallower soil layers, and can 
also be tapped by plant roots directly during dry periods (Fan 
and Miguez-Macho, 2010). Heat exchange beneath the soil 
column affects the seasonal cycle of soil temperature and its 
response time to global warming, which is particularly im- 
portant for modeling permafrost processes (Nicolsky et ah, 
2007; Alexeev et ah, 2007). 

In view of these considerations, a number of studies in re- 
cent years have coupled models of deep water and heat stor- 
ages with models of land- surface and atmosphere processes. 
Many published studies have worked in local to regional do- 
mains and over timescales of days to months. For example, 
Maxwell et ah (2007) showed that including groundwater 
resulted in more realistic spatial soil moisture variability in 
an Oklahoma watershed which, in turn, affected local atmo- 
spheric convection and boundary layer height over the short 
(36 h) study period. Regional model studies by Anyah et ah 
(2008) and by Jiang et ah (2009) both found that ground- 
water storage and flow tends to increase available soil mois- 
ture and hence warm-season evaporation rates in the semiarid 
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central and western United States, with impacts on simulated 
seasonal temperature and precipitation. 

Representations of subsurface water on a global scale have 
mostly been tested in the Community Land Model (CLM) 
developed by the National Center for Atmospheric Research. 
The new versions 3.5 and 4.0 of CLM include such a repre- 
sentation (Oleson et al., 2008; Lawrence et al., 2011), based 
on the SIMGM model of Niu et al. (2007). In SIMGM, an un- 
confined aquifer is parametrized as a single deep layer under- 
lying the model soil layers and exchanging water with them, 
as well as contributing to runoff. Inclusion of this deep layer 
improved the large dry bias in previous versions of CLM 
caused by too-quick drainage and insufficient soil water ca- 
pacity (Bonan and Levis, 2006; Stockli et al., 2008; Pan et al., 
2008; Choi and Liang, 2010). The aquifer in CLM does not 
exchange heat with the shallower soil layers. Instead, there 
is a separate representation of deep bedrock layers for heat 
transfer (Lawrence et al., 2008). This separation between 
heat and water flux is unphysical (since the “deep bedrock” is 
assumed to occupy the same physical space as the “aquifer” 
yet not interact with it), and energy conservation is problem- 
atic since the aquifer temperature is not tracked as a state 
variable. Lo and Famiglietti (2011) showed that the wetter 
soil resulting from the aquifer in CLM 3.5 leads to signifi- 
cant changes in simulated global climate when CLM is cou- 
pled with the Community Atmosphere Model (CAM) 3.5. 
Global land evaporation increased by 9 % and global land 
precipitation by 3 %, with much larger changes over partic- 
ular regions and seasons - for example, over central North 
America summer precipitation increased by some 25 %. Lo 
and Famiglietti (2011) discuss the changes in atmospheric 
water transport that lead to these precipitation responses. 

Here, we report on the impacts of including a simple 
aquifer parametrization in a different atmospheric general 
circulation model (GCM), the Goddard Institute for Space 
Studies (GISS) “ModelE”. Unlike most previous studies, our 
experiments allow us to consider the impacts of aquifer-soil 
heat and water fluxes separately as well as in combination 
in an idealized but physically consistent (mass and energy 
conserving) framework. This study builds on previous work 
investigating soil moisture feedbacks on climate in ModelE 
(Krakauer et al., 2010), which found that soil moisture dy- 
namics affected the persistence and intercorrelation of cli- 
mate variables such as precipitation, evaporation, and sur- 
face air temperature. Our main scientific goal in this paper is 
to examine the impact of aquifer heat and water exchanges 
on the mean climate as well as on the seasonality and inter- 
annual variability of soil temperature and soil moisture and 
of climate variables, such as evaporation, precipitation, and 
surface air temperature which respond to soil conditions. 


2 Methods 

2.1 ModelE 

Our modeling experiments were conducted with the God- 
dard Institute for Space Studies (GISS) ModelE, a state- 
of-the-art atmosphere GCM, run at a resolution of 2° lati- 
tude by 2.5° longitude (Schmidt et al., 2006; Hansen et al., 
2007). The current land surface representation in ModelE 
(Abramopoulos et al., 1988; Rosenzweig and Abramopou- 
los, 1997; Aleinov and Schmidt, 2006) can include separate 
bare soil, vegetated soil, and lake fractions within each grid 
cell. Soil moisture, soil temperature, and land-atmosphere 
fluxes are computed for the two separate soil columns when- 
ever bare soil and vegetated areas coexist within a grid cell. 
Vegetation and phenology are based on Matthews (1983) and 
Matthews (1984), while photosynthesis and stomatal con- 
ductance are computed using the functions of Farquhar et al. 
(1980) and Ball et al. (1987). Each soil column is discretized 
into six soil layers with a total depth of 3.5 m, with runoff 
from each layer depending on its water content and soil prop- 
erties. In the current configuration, there are no heat or water 
fluxes across the bottom of the soil column. 

Soil properties are derived from the data sets described in 
Rosenzweig and Abramopoulos (1997), except that bedrock 
fractions in the deepest soil layers (layers 5 and 6) were re- 
duced to the layer 4 fraction. Vegetation characteristics (leaf 
area index, maximum vegetation height, and root density and 
depth) are also as specified in Rosenzweig and Abramopou- 
los (1997), except for changes to root density and depth. The 
roots of woodland, deciduous tree, and evergreen tree veg- 
etation types are redistributed, so that 20 % of the roots are 
in the bottom two layers (below 1 m). These modifications to 
bedrock function and root depth allow vegetation greater ac- 
cess to deeper soil moisture, and permit more extensive inter- 
action between the soil and the aquifer layer. Multiple stud- 
ies have identified access to deeper soil moisture as an im- 
portant aspect of water uptake by plants during dry periods; 
deep roots have been found to take up substantial amounts of 
water from deep, wet soil when shallow layers dry out, even 
if their contribution to root biomass is small (Jackson et al., 
2000; Feddes et al., 2001; Lee et al., 2005; El Maayar et al., 
2009). 

2.2 Aquifer representation 

We developed a simple, computationally efficient, water- and 
energy-conserving implementation of water and heat flow 
vertically between the soil column and an unconfined aquifer. 
The aquifer is represented as a single additional layer under- 
lying the bottom soil layer. 

The state variables added were the aquifer water content 
w aq (volume of liquid water equivalent per unit area, m) and 
heat content h aq (energy per unit area, J m -2 ). We assumed 
hydraulic and thermal equilibrium within the aquifer layer 
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(i.e., uniform water potential and temperature). The propor- 
tionality factor between change in aquifer water content and 
change in the equilibrium water table is the specific yield S y : 


Az a 


AWaq 

Sy ' 


( 1 ) 


In reality, S y depends on the water table depth za as well 
as on the profile of the specific moisture capacity above the 
water table. We simplified it by assuming a constant S y = 
0.2, which Niu et al. (2007) found to give reasonable results 
globally. With constant S y , we can write 
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( 2 ) 


where ziower is the bottom depth of the aquifer, taken to be 
— 53.5 m. A full aquifer corresponds to za equaling the depth 
of the bottom soil layer, —3.5 m. 

The aquifer heat content h aq is equal to the aquifer water’s 
specific enthalpy 77 a q times its density p times its volume 
w aq , plus the corresponding values for the bedrock (assumed 
to occupy 1 — S y of the aquifer volume and to have a specific 
heat capacity Hp of 2.4 x 10 6 J m -3 K -1 ). The heat content 
yields a corresponding aquifer temperature T aq (cf. Rosen- 
zweig and Abramopoulos, 1997, Eq. 2). 

We assume that at the bottom of the modeled soil lay- 
ers, the water potential (in units of m) is at its value for the 
aquifer, which is za- At the midpoint of the bottom model 
soil layer (layer 6), which has total thickness Azn = 1.53m 
and (vertical) hydraulic conductivity TGv, the water potential 
is equal to the model-calculated value for the bottom layer, 
hu- A finite-difference approximation to Darcy’s law gives 
for the water flux q between the soil and aquifer (volume of 
liquid water per unit area per unit time): 


^soil-aquifer — 


K n 


2(h N - za) 
a ZN 


( 3 ) 


Thus there is no vertical water flow in or out the aquifer if 
z A = /*jv, while there can be flow from the aquifer into the 
soil if za > htf. 

Note that unlike in SIMGM, there is no runoff directly 
from the aquifer. ModelE has subsurface runoff originat- 
ing from the soil layers depending on their water content 
(Abramopoulos et al., 1988). Total runoff remains roughly 
constant after adding the aquifer layer. 

Heat flux q n between the aquifer and soil (W m -2 ) is given 
by 


2(7 a — T a q) 

^H; soil-aquifer = N I - <7 soil-aquifer CpT, (4) 


where K h ; ;v is the bottom soil layer’s (vertical) thermal con- 
ductivity (W m -1 K -1 ), 7V is the temperature of the bottom 
soil layer, and in the second term the water specific heat C p 


and the temperature T are evaluated on the upstream side (cf. 
Rosenzweig and Abramopoulos, 1997, Eq. 1). 

Water and heat flows given by these expressions for each 
timestep are restricted to ensure numerical stability (specifi- 
cally, the aquifer water content is not allowed to decline by 
more than half at each timestep) and to prevent the aquifer 
from overfilling (maximum aquifer water content, given the 
parameter values we chose, is 10 m) or emptying (minimum 
water content set at 10 -6 m). Also, the hydraulic conductiv- 
ity is reduced to zero if the water on the upstream side is 
frozen. 

As previously in ModelE and in other climate mod- 
els (Nicolsky et al., 2007), geothermal heat flow is ne- 
glected in assuming no-flux bottom boundary conditions. 
Adding geothermal heat flow would be expected to warm our 
aquifer layer by up to about 1 K, taking a typical continental 
geotherm of 25Kkm -1 (Mock et al., 1997). This heat flow 
may be more important over multi-millennium simulations, 
particularly for representing conditions at the bottom of ice 
sheets (Pollard et al., 2005). 

2.3 Model runs and analysis 

Two ModelE runs were conducted starting from the same ini- 
tial conditions, one (Control) without an aquifer and with the 
usual no-flux boundary condition at the bottom of the soil 
layers and the other (Aquifer) with an aquifer layer as de- 
scribed. Both runs included dynamic atmosphere and land 
surface processes. All climate forcings were set to 1850 con- 
ditions, and climatological sea surface temperatures and sea 
ice were specified based on 1876-1885 average conditions 
from the HadlSST 1.1 data set (Rayner et al., 2003). For the 
Aquifer run, the aquifer was initialized with heat and water 
contents that corresponded to the monthly mean temperature 
and water potential of the bottom soil layer taken from the 
end of a previous run. 

Two additional ModelE runs, WaterOnly and HeatOnly, 
were carried out in order to evaluate whether particular dif- 
ferences found between the Control and Aquifer runs could 
be attributed primarily to the heat flux or to the water flux 
between the aquifer and the soil column. Both these runs 
started from the same initial condition as the Aquifer run and 
used the same model formulation, except as follows. In Wa- 
terOnly, the thermal conductivity in Eq. (4) was set to zero; 
therefore, there was no heat transfer between the aquifer and 
the soil column (except for heat carried by water fluxes). In 
HeatOnly, the hydraulic conductivity in Eq. (3) was set to 
zero; thus there was no transfer of water between the aquifer 
and the soil column, and the aquifer water content stayed 
fixed at its initial value. 

The Aquifer, HeatOnly, and WaterOnly runs were each run 
to approximate equilibrium. In the Aquifer and WaterOnly 
runs, trends in the aquifer water content were extrapolated 
several times by fitting exponential decay functions to the 
time series and then re-initializing runs with the extrapolated 
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Fig. 1. (a) Mean aquifer water content (expressed as fraction of capacity) in the Aquifer run. (b) Standard deviation (mm) of monthly mean 
aquifer water content in the Aquifer run. 


water contents in an attempt to speed up convergence. Time 
series of hemispheric and global means in model outputs in- 
cluding surface temperature, precipitation, evaporation, soil 
temperature and moisture, and aquifer temperature and water 
content were inspected to verify that all runs were in approx- 
imate steady state. All results shown here are based on the 
last 50 yr from each simulation. 

In the Aquifer and WaterOnly runs, the equilibrium 
aquifer water content was primarily controlled by climate, 
with a full aquifer in moist areas such as the Congo basin and 
western Europe and an empty aquifer in deserts (Fig. la). As 
might be expected, under both these extreme conditions the 
aquifer water content showed little time variability, while the 
greatest variability (indicative of significant aquifer-soil wa- 
ter fluxes) was seen in areas of intermediate dryness such as 
the central United States, the Eurasian steppes, and northeast- 
ern China, where the aquifer water content was intermediate 
and fluctuating (Fig . lb). 

Soil temperature is shown (in the Tables and in Fig. 3) av- 
eraged across the soil column (to 3.5 m depth), and soil water 
content is similarly summed across the soil layers. The sig- 
nificance of differences between Control and the other runs 
(Aquifer, WaterOnly, HeatOnly) was estimated using Stu- 
dent’s t test on time series of the difference between two 


runs taken at annual resolution, with the degrees of free- 
dom adjusted based on the observed lag-1 autocorrelation 
of the time series. Global mean interannual standard devi- 
ation of a climate variable was defined as the average of the 
interannual standard deviation across grid cells, and global 
mean lagged autocorrelation was defined as the autocorre- 
lation of monthly-mean climate variables at a specified lag 
(3, 12, or 48 months), averaged across grid cells. For the 
interannual standard deviations and autocorrelations, stan- 
dard errors of differences in the global mean between runs 
was calculated using their variances and covariances on the 
grid scale as estimated by fitting a spatial Gaussian random 
process with exponential variogram (Cressie, 1993) to the 
map of differences. 

3 Results 

3.1 Aquifer impact on mean climate and seasonal cycles 

As Table 1 shows, soil-aquifer fluxes had remarkably lit- 
tle effect on the mean of climate variables in ModelE. Soil 
moisture increased when soil-aquifer water fluxes were in- 
cluded, because the aquifer provided a storage for excess wa- 
ter when the soil is saturated that can be released as the soil 
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Fig. 2. Month-latitude climatologies of surface air temperature, (a) Control run. (b) Aquifer minus Control run. (c) WaterOnly minus Control 
run. (d) HeatOnly minus Control run. For (b-d), differences not significant at the 0.05 level are hatched. 


HW-CTRL (Soil Temp, K) 


HW-CTRL (Soil Moisture, mm) 



Fig. 3. Month-latitude climatologies, Aquifer minus Control run. (a) Soil temperature, (b) Soil moisture, (c) Evaporation, (d) Precipitation. 
Differences not significant at the 0.05 level are hatched. 
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Table 1. Mean of climate variables 



Control 

Aquifer 

WaterOnly 

HeatOnly 

Soil temperature (°C) 

14.30 

14.29 

14.30 

14.30 

Soil water (mm) 

687.4 

690.6** 

690.5* 

689.1 

Surface air temperature (°C) 

14.12 

14.15 

14.13 

14.16** 

Evaporation (mm day - 1 ) 

1.996 

1.996 

1.995 

1.996 

Precipitation (mm day - 1 ) 

2.982 

2.979 

2.979 

2.986 


Averaged over land grid cells (excluding ice caps). Significance level of differences between Control and 
other runs: *0.05, **0.01. 


Table 2. Interannual standard deviation of climate variables. 



Control 

Aquifer 

WaterOnly 

HeatOnly 

Soil temperature (°C) 

0.49 

0.24** 

0.47* 

0.24** 

Soil water (mm) 

32.5 

26.8** 

27.3** 

31.7* 

Surface air temperature (°C) 

0.56 

0.52* 

0.52* 

0.52* 

Evaporation (mm day - 1 ) 

0.134 

0.132 

0.130 

0.132 

Precipitation (mm day - 1 ) 

0.365 

0.363 

0.358 

0.363 


Averaged over land grid cells (excluding ice caps). Significance level of differences between Control and 
other runs: *0.05, **0.01. 


dries, but this increase was quantitatively tiny (less than 1 %) 
and did not lead to significant change in evaporation or sur- 
face air temperature. Surface air temperature increased when 
soil-aquifer heat fluxes were included, perhaps due to the in- 
teraction of seasonal heat fluxes (see below) with the snow 
thermal rectifier (Goodrich, 1982) at high latitudes, but this 
increase is again tiny (less than 0.1 K). 

Considering changes by season as well as latitude band 
tells a more interesting story. The soil-aquifer heat flux 
greatly dampens the amplitude of the seasonal cycle in tem- 
perature in the deeper soil layers, which lags the seasonal 
cycle in surface insolation by some 2-3 months. This heat 
flux propagates to the surface with additional phase lag, with 
the result that the surface warms in February-August (by up 
to ~ 0.7 K at high latitudes) and cools in September-January 
(by similar amounts) in the Aquifer and HeatOnly simula- 
tions, actually increasing the seasonal amplitude of tempera- 
ture at the surface (Fig. 2). The change in soil temperature is 
greater than in surface air temperature, with warming/cooling 
zonally exceeding 2 K (Fig. 3a). 

Soil moisture considered by latitude shows significant, 
though modest, changes (typically 1-3 %) (Fig. 3b). Notably, 
soil moisture decreases around 50° N, possibly because of 
greater evaporation (Fig. 3c) during warmer springs and sum- 
mers (Fig. 2b). Soil moisture increases further north, as in- 
creasing warm-season evaporation there is more than bal- 
anced by increasing summer precipitation (Fig. 3d) associ- 
ated with warmer conditions. 

The reduced amplitude of the soil temperature seasonal 
cycle means that more grid cells are modeled to have 
permafrost, which we operationally define as the bottom 
model soil layer (layer 6, 1.97-3.50 m depth) never warming 


above freezing during a run. The permafrost area almost 
doubles from 9.1 million to 16.1 million km 2 , in better 
agreement with observation-based inventories which yield 
a range of 12.2-17.0 million km 2 (Zhang et al., 2000). 
Permafrost expands south and east in Siberia, across the 
Tibetan plateau, and around the Arctic coast in North 
America (Fig. 4; compare with the observation-based 
permafrost map at http://www.grida.no/graphicslib/detail/ 
permafrost- extent- in- the- northern-hemisphere _ 1 26 6) . 

Given the importance of permafrost for sequestering carbon, 
this change might be expected to have significant impacts on 
modeled biogeochemistry and response to global warming 
(Koven et al., 2011). 

3.2 Aquifer impact on climate variability and 
persistence 

Table 2 shows the change between runs of the mean inter- 
annual variability in the climate variables previously consid- 
ered. Soil temperature and moisture are both stabilized by 
adding the aquifer layer. As expected, soil temperature re- 
sponds more to the soil-aquifer heat flux, while soil mois- 
ture responds more to the water flux, but both fluxes have 
some effect on the variability in both quantities. Surface air 
temperature variability is also reduced, again showing that 
the influence of aquifer-soil fluxes is detectable in the at- 
mosphere, but evaporation and precipitation do not show 
significant change. 

Table 3 shows the change between runs in the lagged au- 
tocorrelation of climate variables at timescales of 3, 12, and 
48 months, a measure of the climate system’s memory at 
seasonal to interannual timescales. Soil temperature shows 
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Fig. 4. Permafrost extent in the Control run (light blue) and additional permafrost area in the Aquifer run (dark blue). 


much reduced autocorrelation due to soil-aquifer fluxes (pri- 
marily the heat flux) at 3 -month lag, while at 1-year lag there 
is little autocorrelation and no significant difference between 
runs. Soil moisture shows stronger autocorrelation, which 
is decreased especially by the soil-aquifer water flux; the 
aquifer acts as a reservoir with much longer timescale than 
1 yr, and thus dampens the land-surface memory at annual 
and shorter timescales. This loss of memory is also reflected, 
to a lesser but still significant extent, in decreased autocorre- 
lations of surface air temperature (at 3 -month lag) and evap- 
oration (at a lag of 1 yr). At a 4yr lag, all autocorrelations 
are close to zero and few of the differences between runs are 
significant. 

4 Discussion 

We have used an idealized model of soil-aquifer heat and 
water fluxes to assess the likely impact of those fluxes on 
different aspects of global climate. To what extent might our 
conclusions change as our aquifer model becomes more com- 
plex and realistic? 

For heat flux, it is easy to calculate the effect of aquifer 
depth and vertical discretization, assuming that thermal con- 
ductivity is approximately constant with depth. A one-layer 
aquifer will somewhat overstate seasonal soil-aquifer heat 
fluxes compared to a model with finer vertical discretization 
(Alexeev et al., 2007), so that the effect we find on climate of 
the soil-aquifer heat flux can be regarded as an upper limit of 
what we would see as we make our subsurface model more 
realistic. 

Similarly, assuming a uniform water potential within the 
aquifer layer (equivalent to infinite hydraulic conductivity 
within the layer), as we did here, will tend to lead to water 
fluxes that are higher than would be expected in reality for 
the same mean water contents in the aquifer and soil layers. 
Thus the effect we find on climate of the soil-aquifer water 
flux can also be regarded as an upper limit. 

Beyond such discretization effects, which could be mit- 
igated by dividing the aquifer into several thinner layers, 


developing realistic models of subsurface water fluxes is 
complicated by the wide variation in subsurface hydraulic 
properties, which makes it difficult to develop a model that 
gives the right flux magnitudes on the grid scale used in 
GCMs (Choi et al., 2007). Three-dimensional subsurface 
variability in specific yields and hydraulic conductivity, men- 
tioned by Miguez-Macho et al. (2008) as a key source of un- 
certainty in modeling water tables at high horizontal resolu- 
tion, requires subsurface property data sets usable for large- 
scale modeling, similar to those available for surface soil 
properties. An initial step has been taken by Gleeson et al. 
(2011), who compiled representative values for the perme- 
ability of various lithologies based on watershed-scale hy- 
drologic studies. Available water table and streamflow data 
along with total water storage changes from the GRACE 
satellites may also be useful for calibration of parameters in 
GCM simple aquifer representations, using one of several al- 
ready demonstrated approaches (Lo et al., 2010; Werth and 
Guntner, 2010; Becker et al., 2011; Vergnes and Decharme, 
2012), particularly when the GCM runs include more real- 
istic climate variability (compared to the climatological sea 
surface temperatures and greenhouse gas concentrations used 
in the runs reported here, which were intended to represent 
an idealized, stable preindustrial climate). Soil data obtained 
at different depths, now being aggregated in the International 
Global Soil Moisture Network (ISMN, which holds soil tem- 
perature as well as water content data) (Dorigo et al., 201 1), 
may also be useful in validating the modeled effects of wa- 
ter and heat fluxes on the soil temperature and moisture pro- 
files and on their temporal correlation patterns, though mean- 
ingfully comparing global model outputs to in situ measure- 
ments of soil properties remains a challenge (Seneviratne 
et al., 2010; Luo et al., 2012). Borehole temperature mea- 
surements may also be useful in validating model heat fluxes 
over depths up to hundreds of meters under conditions of cli- 
mate change (Beltrami et al., 2006). 

For water flux, horizontal as well as vertical flow within 
the subsurface is likely to be important at least in setting 
equilibrium water distributions, such as relatively high water 
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Table 3. Autocorrelation of climate variables at 3, 12, and 48 month lags. 



Control 

Aquifer 

WaterOnly 

HeatOnly 

3 month lag: 

Soil temperature (°C) 

0.449 

0.151** 

0.428* 

0.145** 

Soil water (mm) 

0.595 

0.578* 

0.574* 

0.592 

Surface air temperature (°C) 

0.055 

0.036* 

0.035* 

0.030* 

Evaporation (mm day - 1 ) 

0.081 

0.077 

0.075 

0.075* 

Precipitation (mm day - 1 ) 

0.010 

0.005 

0.003* 

0.004 

12 month lag: 

Soil temperature (°C) 

0.053 

0.029 

0.042 

0.045 

Soil water (mm) 

0.272 

0.229** 

0.235** 

0.258* 

Surface air temperature (°C) 

-0.003 

-0.018 

-0.014 

-0.004 

Evaporation (mm day - 1 ) 

0.010 

0.000** 

-0.001** 

0.002* 

Precipitation (mm day - 1 ) 

-0.016 

-0.021 

-0.019 

-0.014 

48 month lag: 

Soil temperature (°C) 

-0.011 

0.007 

-0.024 

0.015 

Soil water (mm) 

0.071 

0.055 

0.048** 

0.068 

Surface air temperature (°C) 

-0.026 

-0.012 

-0.018 

-0.011 

Evaporation (mm day - 1 ) 

-0.017 

-0.019 

-0.013 

-0.009* 

Precipitation (mm day - 1 ) 

-0.021 

-0.022 

-0.019 

-0.015 


Averaged over land grid cells (excluding ice caps). Significance level of differences between Control and 
other runs: *0.05, **0.01. 


tables in valleys, and also for sustaining baseflow in rivers 
(Toth, 1963; Fan et al., 2007). Horizontal subsurface flow 
may also be important as a driver of long-timescale hydrocli- 
mate persistence (Bierkens and van den Hurk, 2007). Again, 
representing this at the GCM grid scale is an unsolved prob- 
lem, and will probably require estimating effective horizon- 
tal hydraulic conductivities and rates of subsurface runoff 
from groundwater simulations run and validated at much 
higher resolution (Fan et al., 2007, 2013; Lam et al., 2011; 
Sutanudjaja et al., 201 1). 

Another way in which our model could be refined is to 
include a fraction of deep roots that reach below the 3.5 m 
deep soil column (Stone and Kalisz, 1991), whether in per- 
meable deep soil or regolith (Nepstad et al., 1994; Markewitz 
et al., 2010) or even in fractured bedrock (Schwinning, 2010; 
Roering et al., 2010). This direct connection of transpiring 
vegetation to “aquifer” layers would likely increase their im- 
pact on water budget variables such as evaporation (Fan and 
Miguez-Macho, 2010). 

Our results show that soil-aquifer interactions affect the 
seasonality and persistence of climate variables as well as, 
to a small extent, their mean (e.g., soil moisture increases by 
0.5 %). However, the impacts that we find of aquifer inter- 
actions on mean climate are much smaller than those shown 
by Lo and Famiglietti (2011) for CLM, where, for example, 
including the aquifer component led to increases of 9 % in 
global mean land evaporation and 3 % in global land pre- 
cipitation, while we see changes of much less than 1 % in 


both these variables. Plausibly, this is because adding the 
aquifer parametrization to CLM dispensed with the prior 
free-drainage bottom boundary condition for soil moisture, 
which had led to chronically dry soils (Zeng and Decker, 
2009). ModelE, by contrast, has a no-flow bottom bound- 
ary condition in the Control configuration, which leads to 
the same result (no water flux at the soil lower boundary) 
as an aquifer at hydraulic equilibrium with the deep soil. 
Thus adding the aquifer modified the variability of soil tem- 
perature and moisture in ModelE but had much less impact 
on the mean climate state than in CLM. Based on our re- 
sults, we hypothesize that modifying a land surface model 
with a no-flux bottom boundary condition to include soil- 
aquifer heat and water fluxes that are physically consistent 
and of realistic magnitude will in general lead to very lim- 
ited change in mean temperatures and hydrological budgets, 
although the changes in seasonal and interannual variability 
might be quite significant from, for example, a biogeochem- 
ical standpoint (as for permafrost extent). 

By controlling heat and water fluxes separately, we found 
that both can impact variability in such quantities as soil 
moisture and evaporation. Along similar lines, soil moisture 
changes can significantly influence soil temperature profiles 
and permafrost extent by changing the effective thermal con- 
ductivity (Subin et al., 2013). These results underline the de- 
sirability of simulating soil-aquifer heat and water fluxes to- 
gether in a physically consistent way, since one or the other 
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in isolation will not fully represent the impact of the deep 
subsurface on climate variability and persistence. 

We found that soil-aquifer fluxes tend to decrease the au- 
tocorrelation of land-surface and climate variables (soil tem- 
perature and moisture, surface air temperature, evaporation, 
and precipitation) on timescales of up to 1 yr. On longer 
timescales, it is plausible that aquifers contribute to multiyear 
climate memory and persistence (e.g., long droughts and plu- 
vials), although we did not find conclusive evidence for this 
in the analyses presented here. One caveat for interpreting 
our model results is that our runs had artificially low interan- 
nual variability in climate forcing and sea- surface boundary 
conditions (SSTs kept at climatology, constant greenhouse 
gas levels); a next step will be to examine how climate persis- 
tence changes in the presence of more variability, as well as 
in the presence of changing climate forcings such as green- 
house gas levels. As another example of long-term change, 
drawdown of groundwater for irrigation, now observable in 
many parts of the world (Tiwari et al., 2009; Kustu et al., 
2010; Famiglietti et al., 2011), may contribute to regional 
drying, although the added irrigation results in moister soil, 
higher evaporation, surface cooling, and higher precipitation 
over the short term (Sacks et al., 2009; Cook et al., 2011; 
Ferguson and Maxwell, 2012). 

5 Conclusions 

We present and test an idealized parametrization of water 
and heat flow between a subsoil aquifer and the soil column 
in the GISS ModelE GCM. The simple aquifer parametriza- 
tion used here suffices to show the impact of these fluxes on 
seasonality and the potential impact of groundwater on per- 
sistence of soil moisture and climate. Many aspects of this 
parametrization (e.g., single aquifer layer, no flow between 
grid cells, no roots below the soil column) can usefully be 
refined to study these impacts with more detail and accuracy, 
especially as they relate to temperature and moisture changes 
resulting from anthropogenic climate change and aquifer de- 
pletion. 
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